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EIGENMATRIX  REPRESENTATIONS  OF  RADIANCE  DISTRIBUTIONS 
IN  LAYERED  NATURAL  WATERS  WITH  WIND-ROUGHENED  SURFACES 


Rudolph  W.  Preisendorfer 


ABSTRACT.  This  report  develops  analytic,  closed-form  solutions 
for  radiance  distributions  in  natural  waters  such  as  lakes  and 
seas.  The  solutions  are  valid  in  layered  water  bodies  for  which 
each  layer  has  inherent  optical  properties  (absorption  and 
scattering  functions)  which  are  independent  of  depth  within  that 
layer.  The  water  body  is  assumed  free  of  internal  light 
sources.  The  effects  of  a  wind-blown  air^sea  surface  are 
included.  This  work  extends  to  the  radiance  level  certain  results 
which  were  previously  known  to  hold  for  irradiances.  The 
eigenmatrix  formalism  developed  here  is  convenient  for  numerical 
computation  of  radiance  distributions,  given  the  inherent  optical 
properties  of  the  water  and  the  desired  boundary  conditions  at  the 
water  surface  and  bottom  (the  direct  problem).  Moreover,  the 
formalism  suggests  an  algorithm  for  solving  the  inverse  problem: 
the  determination  of  the  inherent  optical  properties  from 
measurements  of  the  radiance  distribution  within  a  water  body. 


1.  INTRODUCTION 

We  develop  here  a  method  for  solving  the  equation  of  transfer  for 
radiance  in  a  piecewise  homogeneous,  source-free,  plane  parallel  water  body 
with  a  wind-roughened  air-water  surface.  The  assumption  of  homogeneity  within 
each  layer  of  the  medium  allows  a  closed-form  type  of  solution  which  views  the 
radiance  distribution  at  each  depth  within  a  layer  as  a  linear  combination  of 
purely  exponential  elementary  components  which  decrease  or  increase  with 
depth.  This  mode  of  decomposition  of  the  multi-directional  light  field  has 
the  same  simple  visualizability  as  the  classic  two-flow  irradiance  model  of 
the  light  field  from  which  the  theory  of  radiative  transfer  in  stratified 
media  began  (cf.  Schuster,  1905).  Moreover  the  new  representation  allows 
explicit  analytic  and  algebraic  formulas  to  be  developed  for  such  basic 
properties  of  the  medium  as  the  radiance  reflectance  and  radiance 
transmittance  of  the  component  layers,  the  radiance  distribution  within  each 
layer,  and  also  the  asymptotic  radiance  distribution  evolving  in  an  infinitely 


deep  homogeneous  lower  layer  of  a  medium.  Numerical  solutions  are  readily 
forthcoming  from  such  algebraic  and  analytic  representations,  and  the  boundary 
conditions  for  a  wind-ruffled  air-water  surface  endow  the  model  with  an 
ability  to  handle  realistic  reflectance  and  transmittance  activity  at  the 
surface. 

The  present  study  builds  on  an  earlier  work,  Mobley  and  Preisendorfer 
(1988),  which  describes  in  detail  the  so-called  Natural  Hydrosol  Model  (NHM) 
for  the  computation  of  radiance  distributions  in  natural  waters  with  wind¬ 
blown  surfaces.  (The  NHM  does  not  require  the  assumption  of  piecewise 
homogeneity  of  the  water  body.)  In  the  interests  of  brevity,  we  shall  draw 
freely  from  the  results  of  this  previous  work;  familiarity  with  Mobley  and 
Preisendorfer  (1988)  is  therefore  a  prerequisite  for  the  full  understanding  of 
the  present  work.  The  essential  feature  of  the  NHM  is  the  quad-averaging  of 
the  radiance  distribution  over  quadrilateral  subsets  of  the  unit  sphere  of 
directions,  and  the  exact  splitting  of  the  azimuthal  structure  of  the  quad- 
averaged  radiance  distribution  into  spectral  modes  by  harmonic  analysis.  This 
double  decomposition  of  the  radiance  distribution  has  the  effect  of  reducing 
the  integrodif f erential  equation  of  transfer  for  radiance  to  a  set  of  coupled 
ordinary  differential  equations,  one  to  each  azimuthal  spectral  mode.  From 
this  point  on,  the  techniques  of  the  linear  interaction  principle  (cf. 
Preisendorfer,  1976,  vol.  IV)  can  be  applied  to  the  family  of  differential 
equations  associated  to  each  spectral  mode. 

In  homogeneous  media,  the  system  matrices  of  each  spectral  family  of 
differential  equations  have  depth-independent  entries;  this  allows  an 
extremely  useful  algebraic  decomposition  of  tne  system  matrix  into  its 
eigenvalues  and  eigenvectors  in  each  homogeneous  layer  of  the  medium.  The 
eigenvalues  of  the  system  matrix  turn  out  to  be  the  desired  exponential  modes 


of  decay,  and  the  eigenvectors  become  the  framework,  for  the  directional 
structure  of  the  radiance  distributions. 

For  those  who  wish  to  see  the  ground  from  which  the  present  eigenmatrix 
procedure  has  sprung,  we  include  in  Appendix  A  the  differential  equations  of 
the  two-flow  model  of  the  irradiance  field,  the  modern  descendants  of 
Schuster's  (1905)  equations.  It  is  indicated  how  the  classic  two-flow  model 
establishes  an  algebraic  pattern  that  generalizes  from  the  simple  upward  and 
downward  decomposition  of  photon  flows  to  the  multitude  of  photon  flows  in  the 
present  context.  For  those  who  are  coming  upon  the  notions  of  the  linear 
interaction  principle  and  its  associated  invariant  imbedding  procedures  for 
the  first  time,  the  exposition  in  Appendix  A  should  perhaps  be  studied  before 
going  on  to  section  2,  below. 

Sections  2  through  5  merely  restate  various  results  which  are  rigorously 
developed  in  Mobley  and  Preisendorfer  (1988).  The  goal  of  this  review  is  the 
local  interaction  principles  as  expressed  in  equation  (5.16),  along  with  the 
associated  boundary  conditions  (5.18),  (5.19)  and  (5.22).  The  real  work,  of 
this  report  then  begins  in  §6.  A  high  point  (at  least  for  the  author)  in  the 
present  exposition  comes  in  sections  7  and  8,  where  the  physical  meanings  of 
the  eigenvalues  of  the  system  matrices  become  clear,  and  where  the 
eigenvectors  of  the  system  matrices  are  shown  to  give  rise  to  linearly 
independent  pieces  of  the  radiance  distribution.  The  development  continues 
until  the  main  goals  of  this  work  are  reached  in  §16  and  §17. 

Acknowledgments .  Initial  numerical  explorations  of  the  eigenvalue  and 
eigenvector  problems  defined  here  were  made  by  Curtis  D.  Mobley,  who  was 
partially  supported  by  a  TOGA  (Tropical  Ocean  Global  Atmosphere)  Council 
contract,  and  who  also  was  partially  supported  by  the  National  Climate  Program 


Figure  1  shows  the  geometric  arrangement  of  the  body  X[x,z],  x  <  z,  of 


the  natural  hydrosol,  its  upper  boundary  X[a,x],  and  its  lower  boundary 
X[z,b].  In  the  present  discussion  X[a,x]  is  the  infinitesimally  thin  average 
plane  of  the  random  air-water  surface.  We  imagine  level  a  to  be  just  above 
and  level  x  to  be  just  below  the  surface.  X[z,b]  can  be  either  an 
infinitesimally  thin,  opaque,  matte  reflecting  bottom  at  a  finite  optical 
depth  z  below  x,  or  an  optically  infinitely  deep,  homogeneous  medium  below 
depth  z,  with  b  =  ®.  Both  cases  will  be  considered  below. 

The  boundary  conditions  for  (2.1)  at  the  random  upper  air-water  surface 

are 


N(a',£)  =  J  N(x;£')  t(x,a‘,£';£)  dfl(£'  ) 

+  J  N(a;£')  r(a,xj£';£)  dfl(£' )  ,  5  e  E+  (2.4) 

N(x;£)  =  J  N(a-,£' )  t(a,x;c';£)  dft(£' ) 

+  J  N(x;£')  r(x,a;£’;£)  dQ(e')  ,  5  e  =_  (2.5) 


Here  E+  and  H_  are  the  upper  and  lower  hemispheres  of  H,  respectively. 

The  r  and  t  functions  describe  how  radiance  is  on  average  reflected  and 
transmitted  by  the  boundary  surface.  These  are  determined  by  a  Monte  Carlo 
method  as  developed  in  Preisendor f er  and  Mobley  (1985,  1986)  and  applied  to 
the  Natural  Hydrosol  Model  in  Mobley  and  Preisendorf er  (1988). 

The  boundary  condition  for  (2.1)  at  the  bottom  level  z  is  given  by 

N(z*,£)  =  J  N(z-,£' )  r(z,b;£' ;£)  dn(£*  )  ,  £  c  (2.6) 
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3.  QUAD-AVERAGED  EQUATION  OF  TRANSFER  AND  BOUNDARY  CONDITIONS 

The  unit  sphere  5  is  shown  in  Fig.  2  partitioned  into  a  set  of  quads 

which  includes  as  special  cases  a  pair  of  polar  caps.  There  are  m  quad  zones 

above  and  m  quad  zones  below  the  equator,  with  each  cap  counting  as  a  special, 

zone.  Each  hemisphere  is  divided  into  2n  azimuthal  sectors.  In  Fig.  2,  m  *  5 

and  2n  *  20. 

A  non-poLar  or  regular  quad  Quv  is  indexed  by  a  pair  of  integers  u,v 
where  u  =  i,...,m-l  is  the  zonal  index  and  v  =  l,...,2n  is  the  azimuthal 
index.  Regular  quads  Quv  have  equal  angular  widths  Aq> v  =  A$  =  n/n  and 
arbitrary  heights  dpu.  Quad  Quv  is  centered  at  azimuth  angle  $v  =  (v-l)A4>  and 
subtends  a  solid  angle  Quv  ■  A$Auu<  Polar  caps  are  denoted  by  "Qm"  and 
subtend  solid  angles  of  size  =  2itAum.  It  will  be  clear  from  the  special 
notation  developed  below  which  hemisphere  (=+  or  =_)  a  cap  or  quad  is  in. 

The  quad-averaged  radiance  over  quad  Quv  is  defined  by  writing 


N(y;u,v)  =  fl-i  J  N(y;£)  dflU) 


(3.1) 


=  Q_  1  f  J"  N(y;y,<>)  dpd$ 


Here  we  used  an  alternate  description  of  5  via  its  zenith  (u  =  cos9)  and 
azimuth  ($)  coordinates.  When  Quv  is  in  =+  or  in  S_  we  will  write  the 
associated  quad-averaged  radiance  as  "N+(y;u,v)"  or  "N-(y;u,v)", 
respectively.  In  each  of  these  cases,  u  runs  from  1  to  m-1  and  v  from  1  to  2n 
for  regular  quads.  The  quad-averaged  radiance  over  a  polar  cap,  for  which 
u  =  m  and  v  is  undefined,  is  denoted  by  "N+(y ;m, • )"  or  "N-(y;m, • )". 

The  quad-averaged  phase  function  is  defined  by  writing 


Figure  2. — An  example  partitioning  of  the  unit  sphere  into  quads  and  polar 

caps,  for  the  case  of  m  =  5  and  n  =  10.  The  coordinate  system  and  quad 
indexing  scheme  is  centered  in  the  wind  direction.  Quad  Qrg  is  shown  in 
the  lower  hemisphere  of  directions,  E_,  and  Quv  is  shown  in  the  upper 
hemisphere,  H+. 
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N(a;u,v)  =  £  £  N(x;r,s)  t (x,a}r, s  | u, v)  •*■  £  £  N(a;r,s)  r(a,x}r,s |u,v) 


Q  in  = . 

uv  + 


(3.5) 


N(x;u,v)  =  £  £  N(a;r,s)  t(a,x;r,s ju,v)  +  £  £  N(x;r,s)  r(x,a;r,s |u,v) 


Q  in  E 
uv  - 


(3.6) 


where  the  four  surface  transfer  functions  t(a,x),  r(x,a),  t(x,a),  and  r(a,x) 
are  all  defined  following  the  general  pattern 


f(r,s|u,v)  =  Q_1  J  J*  dyd<t>  J  J*  du ’ d<t> '  f  (y '  jy,4>) 

UV  Q  Q 

_  _  uv  r  s 

<W  ^uv  in  = 


(3.7) 


where  "f (y' , <t> *  ;y,<t>)"  denotes  any  of  the  four  transfer  functions  in  (2.4)  or 
(2.5).  The  bottom  boundary  condition  in  quad-averaged  form  is 


N(z;u,v)  =  £  £  N(zjr  ,s)  r(z,b;r,s|u,v) 
r  s 

Quv  in 


(3.8) 
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4.  QUAD- AVERAGED  LOCAL  INTERACTION  PRINCIPLES 

The  quad-averaged  equation  of  transfer  (3.4)  may  be  split  into  two 

statements,  one  for  the  upward  radiances  N+(y;u,v)  and  one  for  the  downward 

radiances  N-(y;u,v)  at  each  level  y,  x  <  y  <  z.  The  isotropy  of  the  volume 

scattering  functions  o(y;£';£),  namely  the  property  that  its  values  depend 

only  on  £*•£  and  not  £'  and  £  separately,  considerably  simplifies  the 

structure  of  the  transfer  equation.  Thus  we  can  write 

Q  in  E  and  Q  in  5 
rs  +  uv  + 


"p+(y;r,s |u,v)"  for  p(y;r,s|u,v)  if 


or 


"p  (y?r,s|u,v)"  for  p(y;r,s|u,v)  if 


Q  in  5  and  Q  in  H 
r  s  -  uv 


Q  in  H  and  Q  in 
rs  +  uv 


(4.1) 


or 


Q  in  E  and  Q  in  E . 
rs  -  'uv  + 

Hence  p+  and  p-  respectively  act  like  local  transmittance  and  reflectance 

functions  in  the  body  of  the  medium,  relative  to  the  horizontal  plane  of  the 

equator  of  E.  We  then  find  from  (3.4)  that 


*Pu  ^  N~(y;u,v)  =  -N  (y;u,v)  ♦  w(y)  £  ^  N+(y;r,s)  p~(yjr,s |u,v) 


r  s 


(4.2) 


+  ui(y)  l  l  N  (y;r,s)  p+(y;r,s  |u,v) 


r  s 


u  =  1 , . . . ,m;  v  =  l,...,2n. 


This  is  a  coupled  pair  of  differential  equation  systems.  The  upward 
system  is  obtained  by  taking  all  upper  signs  together.  This  system  describes 
the  evolution  with  optical  depth  y  of  the  upward  radiances  N+(y;u,v).  The 
downward  system  describes  N_(y;u,v).  The  complete  system  (4.2)  constitutes 


the  local  interaction  principles  or  the  local  forms  of  the  principles  of 


invariance.  See  Preisendorf er  (1965,  p.  103),  H.O.,  Vol.  Ill,  p.  4,  and 


Vol.  II,  p.  295.  The  boundary  conditions  (3.5),  (3.6),  and  (3.8)  hold  also 


for  (4.2). 
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5.  SPECTRAL  FORM  OF  THE  QUAD-AVERAGED  LOCAL  INTERACTION  PRINCIPLES 
We  next  split  the  equation  set  (4.2)  into  smaller  groupings  of  dependent 
variables  by  means  of  Fourier  polynomial  analysis.  This  entails  no  loss  of 
information  of  the  radiance  field  but  considerably  facilitates  the  numerical 
solution  of  the  set  (4.2). 

For  fixed  y  and  u,  N~(y;u,v)  is  a  function  defined  on  the  finite  set 
consisting  of  an  even  number  of  integers  v  =  l,...,2n,  corresponding  to 
azimuth  angles  4>v  *  (v-l)Ad,  as  defined,  above.  We  may  then  represent  this 
v-dependence  of  N^yjUjv)  by 


N“(y;u,v)  =  £  [ A” ( y ; u ,  1 )  cos(l<t>  )  ♦  A~(y;u,l)  sin(a$  )] 

1=0  v  v 


(5.1) 


where 


+  ***  +. 

A~(y;u;l)  =  e~l  £  N“(y;u,v)  cos(8.$  ) 

1  v=i  v 

i  —  0  , . . . ,  n 


(5.2) 


+  411  + 

A~(y;u;i)  =  y~ 1  £  N~(y;u,v)  sin(&$  ) 

*  v=l 

l  —  1,..., n- 1 


(5.3) 


These  Ap(y;u,l),  p  =  1  or  2,  are  the  radiance  amplitudes  for  (harmonic)  mode 
t.  The  factors  e,  and  y.  are  given  by 


ijl 

hS 

Si 


I 


2n  if  1=0  or  l=n 

n  if  l~l f • • • (H“l 
0  if  1=0  or  l=n 


(5.4) 


(5.5) 


(  n  if 

Observe  that  y^  =  ea  *  n  over  the  range  l  = 

Since  $v  =  (v-l)ir/n,  we  see  that  sin(l<t>v)  =  0  for  1  =  0  and  1  =  n.  Hence 
we  will  define  Af(y;u;0)  =  0  and  A^yiujn)  =  0  for  u  =  l,...,m.  Also  note 
that 


A^(y;m;0)  =  N  (y;m, • ) 


Aj (y ;m; 1 )  s  0  )  l~l | . . . f n 


(5.6) 


and  that 


A2(y;m;i)  =  0  ,  a=0,...,n 


(5.7) 


•a 


Therefore  all  radiance  amplitudes  for  the  polar  caps  are  zero,  except  for  the 
cosine  amplitude  for  the  zero  azimuthal  index,  1=0.  Holding  y,  r,  s,  and  u 
fixed,  we  can  represent  the  phase  function,  as  a  function  of  v,  in  Fourier 
polynomial  form: 


+  .  r  ± 

p  (y ; r , s | u,v)  =  I  p  (y;r,u;l)  cosl($  ) 
1=0  3  v 

x  <  y  <  2 

r )  u  “  1 9  •  •  •  9  m  I  s  |  v  ■*  1 , .  • .  ,  2n 


1 

id 


where 
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?  A-(y;a) 


=  A j ( y *  1 )  T(y;0  ♦  A2(y;£)  £(y*,£) 


(5.14) 


where 


A2(y;£)  =  [A2(y; 1;£) , . . . ,A2(y;m-l ;a) ] 


(5.15) 


1  * . . .  *  n  1 


The  matrices  p(y;£)  and  x(y;£)  are  either  mxm  or  (m-l)x(m-l),  as  can  be 
inferred  in  each  case  from  the  number  of  components  of  the  A~(y;£)  vectors, 
p  =  1*2.  These  matrices  are  fully  defined  in  Mobley  and  Preisendorfer 
(1988).  What  should  be  noted  here  is  that  x(y;£)  is  a  local  transmittance 
matrix  in  the  sense  that  it  propagates  Ap(y;£)  into  A±(y;l)  over  an 
infinitesimal  increment  Ay  of  optical  depth.  p(y;£)  acts  as  a  local 
reflectance  matrix  in  the  sense  that  it  respectively  converts  A*(y;£)  into 
Ap(yjl)  over  Ay.  In  this  way  the  rising  and  descending  streams  of  photons  in 
an  infinitesimal  layer  X[y,y+Ay]  of  a  medium  X(x*z]  feedback  to  each  other  and 
generate  the  multiply  scattered  radiance  field. 

The  preceding  three  autonomous  sets  of  coupled  differential  equations  all 
fall  into  the  following  general  pattern 


?  A  <y;*>  =  A  (y*fc)  T (y 5 a )  +  A*(y;£)  S(y*£) 
dy  -p  -p  -  -p~ 

x  <  y  <  z 

p  "  1*2  *  £  ~  0 * . . . ,n 


A  (y;£)  =  (A"(y;l;£),...,A  (y;q;£)] 

-p  p  p 

q  =  m-1  or  m 


Thus  the  amplitude  vectors  A”(y;a)  in  (5.16)  are  lxq  and  the  matrices 
t(y;i)  and  p(y;a)  are  qxq,  where  q  is  either  m-1  or  m,  as  the  case  may  be, 
i.e.,  depending  on  which  of  (5.10),  (5.12),  or  (5.14)  we  are  considering. 

Once  the  three  sets  (5.10),  (5.12),  and  (5.14)  are  solved  we  will  have  at 
each  depth  y  exactly  enough  amplitudes  to  construct  the  radiances  N~(y;u,v), 
u  =  l,...,m;  v  *  l,...,2n  via  (5.1).  We  will  have  for  each  flow  (±),  n+1 
cosine  amplitude  vectors  Aj(y;a),  a  =  0,...,n  and  n-1  sine  amplitude  vectors 
A|(y;l) ,  1  *  l,...,n-l,  each  with  m-1  or  m  components,  as  needed. 

The  surface  boundary  conditions  at  X[a,x]  that  go  with  (5.16)  are 


(5.18) 


A+(a;i)  = 

n 

l 

A+(x;k) 

n 

t  (x,a;k|a)  ♦  Y 

A  (ajk) 

r  (a,x;k 

a) 

-P 

k=0 

-p 

p  k=0 

~P 

“P 

A  (x;t)  = 

n 

l 

A  (x;k) 

n 

t  (a,x;k| a)  +  Y 

A+(x;k) 

r  (x,a;k 

a) 

_P 

k*0 

~p 

P  k=0 

~P 

“P 

p  *  1  or 

2  , 

1*0, 

. . .  ,n,  and  k+a  even. 

(5.19) 


These  are  obtained  by  using  the  Fourier  representations  of  each  member  of 
(3.5)  and  (3.6)  and  reducing  to  the  indicated  forms  in  (5.18)  and  (5.19).  The 
azimuthal  (^-behavior )  symmetries  of  the  random  wind-ruffled  air  water  surface 
are  those  of  an  ellipse  (cf.  Mobley  and  Prei sendorf er ,  1988),  which  among 
other  things  require  t  and  £p  to  vanish  when  k+a  is  odd.  Hence  the  sums  in 
(5.18)  and  (5.19)  may  be  restricted  to  those  values  of  k  and  l  for  which  k+a 
is  even. 

The  entries  of  the  four  transfer  matrices  fp,  £  in  (5.18)  and  (5.19)  are 
given  in  Mobley  and  Preisendorfer  (1988).  Observe  that  these  are  mxm 
matrices,  some  of  which  have  zeros  in  their  mth  rows  or  mth  columns  (see 
Tables  1  and  2,  Mobley  and  Preisendorfer,  1988).  The  essential  point  to  note 
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-behavior  of  r(z,b;r,s |u,v) 


)  COSl(lJ)  ) 
s  v 


2n 

l  r(z,b;r,s| 


u,v)  cos(!l<(>  ) 


#2 


xvxvkv  xv  xvxv  yvxyxv  xvw  v\,  yvvvirv*-.  ^ w 
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6.  THE  FUNDAMENTAL  MATRIX  FOR  RADIANCE  AMPLITUDES  IN  HOMOGENEOUS  LAYERS 
We  may  now  proceed  Co  the  main  interest  in  this  study,  the  solution  of 
the  equation  set  (5.16)  when  w(y)  and  p(yJC'»5)  in  (3.4)  and  hence  p(y;l)  and 
i(y?0  are  independent  of  optical  depth  y  in  an  arbitrary  layer  X[x,z]  of  a 
natural  hydrosol.  The  solution  procedure  we  develop  is  independent  of  whether 
p  =  1  or  2  or  what  mode  index  1  =  0,...,n  is  of  current  interest. 

According!  we  can  until  further  notice  drop  both  "p"  and  "i"  from  the 
notation  in  (5.16).  (Both  p  and  l  will  have  to  be  reinstated  in  §12,  for 
example,  when  boundary  conditions  (5.18)  and  (5.19)  are  to  be  used,  and  also 
when  the  final  Fourier  synthesis  of  N~(y;u,v)  in  (5.1)  is  made). 


A.  -tacic  Local  Interaction  Principles 

*/e  will,  in  accordance  with  the  preceding  notational  convention,  now  work, 
with  the  following  streamlined  versions  of  (5.16)  and  (5.17): 


d  + 


A  (y)  *  A  (y)  _t  +  A+(y)  £ 


A  (y)  =  [A~(y}l ) , . . . ,A  ( y ; q ) ; 
x  <  y  <  z 


(6.1) 


<  .2) 


Thus  6  and  x  are  qxq  matrices  with  constant  entries  and  A~(y)  at  each  depth  y 
are  lxq  vectors.  We  may  further  simplify  the  system  (6.1)  and  (6.2)  by 
writing 


A(y)  =  [A  (y),  A  (y)] 


( lx2q) 


(6.3) 


WMlm® 


WWK8B 


K  i 


(2qx2q) 


(6.4) 


so  that  (6.1)  and  (6.2)  become  the  following  version  of  the  local  interaction 
principlesl 


^  A(y)  *  A(y)  K 
x  <  y  <  z 


(6.5) 


B.  The  Constructive  Definition  of  w(x,y) 

Equation  (6.5)  can  be  integrated  at  once  either  numerically  or 
formally.  We  shall  concentrate  here  on  the  former.  Numerically,  one  would 
choose  an  initial  lx2q  vector  A(x)  and  then  march  (6.5)  down  from  level  x  to 
any  level  y,  x  <  y  <  z,  in  the  given  medium  X[x,z].  There  are  many  such 
initial  vectors  A(x)  from  which  one  may  start.  For  example,  one  may  have 
measured  the  radiances  N~(x;u,v),  u  =  l,...,m;  v  =  l,...,2n  at  level  x  just 
below  the  surface.  From  these  radiances  one  can  find  A~(x)  as  shown  in  (5.2) 
and  (5.3);  whence  A(x).  Then  the  amplitude  A(y)  is  found  by  integrating 
(6.5),  with  A(x)  as  initial  vector,  down  to  any  desired  depth  y,  x  <  y  <  z. 

There  is  one  important  set  of  initial  vectors  A(x),  however,  that  leads 
to  the  general  solution  of  (6.5).  This  is  the  set  of  2q,  lx2q  initial  vectors 
stacked  vertically  to  form  the  2qx2q  matrix: 


10  0  •  •  •  0  0 

0  10  •  •  •  0  0 


w(x,x) 


0  0  0 
0  0  0 


1  0 
0  1 


(6.6) 
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I  jjs  ij*  <V', 


li: 


where  I^q  i-3  a  2qx2q  identity  matrix.  Thus  the  initial  lx2q  vectors  form  the 
2q  rows  of  the  matrix  #(x,x). 

Let  mj(y)  be  the  lx2q  vector  soLution  of  (6.5)  for  the  initial  lx2q 
vector  mj(x)  =  [ 0 , . . . ,  1 , . . . ,0 ]  where  all  components  are  zero  except  that  in 
place  j,  1  <  j  <  2q.  Thus  mj(y)  satisfies  ( 6 . 5 ) i 


—  m^(y)  =  m^£  ,  j  »  l,...,2q 


Write 


*(x, y) 


®i(y) 


rn^Cy) 

Vi(y) 


2?o(y) 


(2qx2q) 


(6.7) 


Then  clearly 


Mix, y)  =  M(x,y)  K 


(6.8) 


«(x,x)  =  I0 
-  -2q 


(6.9) 


The  2qx2q  matrix  Mix, y)  in  (6.7)  is  the  fundamental  matrix  of  the  system  of 
differential  equations  (6.5).  The  importance  of  this  system  rests  in  the  fact 
that  if  A(x)  is  any  lx2q  vector,  then  the  depth  dependent  lx2q  vector 
A(y)  =  A(x)  Mix ,y)  is  the  solution  of  (6.5)  with  initial  vaLue  A(x).  This 


MM* 


m 

m 


r  V  ■* 

.-v>\ 


'v-Vfr/ 

v 


'W\ j 
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follows  formally  on  multiplying  (6.8)  from  the  left  by  A(x)  and  reducing  the 
result  to  (6.5).  The  mathematical  reason  for  this  remarkable  property  is  that 


the  2q  vectors  mj(y)  defined  above  are  linearly  independent  for  each  y, 
x  <  y  <  z.  Therefore  any  solution  m(y)  of  (6.5)  at  some  y  is  a  vector  in  the 
space  spanned  by  the  mj(y)  (cf.  Coddington  and  Levinson,  1955,  pp.  68,  69). 
Thus  we  have  the  mapping  property  of  «(x,y): 

A(y)  =  A(x)  M(x, y)  (6.10) 

x  <  y  <  z 

This  mapping  property  also  can  be  written  down  for  the  fundamental 
solution  m( y,z)  of  (6.5)  for  which  M( y,y)  =  ^q  and  where  w(y,z)  is  obtained 
by  integrating  (6.8)  starting  from  level  y.  The  associated  mapping  property 
is  then  A(y)  m( y,z)  =  A(z).  Combining  this  with  (6.10),  noting  that  we  have 
also  A(z)  *  A(x)  M(x,z),  and  letting  A(x)  be  arbitrary,  we  obtain  the  group 
property  of  M(x,z): 

M(x,z)  =  «(x,y)  M(y,z)  (6.11) 

x  <  y  <  z. 

Setting  x  =  z  in  (6.11),  the  inverse  of  «(x,y)  is  found  to  be 

M(y,x)  =  M“l(x,y)  (6.12) 

x  <  y  <  z 

One  may  for  example  determine  «(y,x)  by  integration,  on  starting  at  level  y 
and  integrating  (6.5)  upward  from  y  to  x.  This  would  result  in  an  associated 
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mapping  property  for  «( y,x)  with  initial  amplitude  A(y).  Alternately,  one  may 
invert  M(x,y)  by  algebraic  or  numerical  means  to  find  M(y,x).  Thus  M_1(x,y) 
has  the  physical  interpretation  of  the  mapping  operator  M( y,x),  i.e.,  from 
level  y  to  level  x  <  y.  trora  our  observations  above  on  the  vector  space 
aspects  of  M(x,y)  it  is  clear  that  M~l(.x,y)  exists  in  all  natural  hydrosols 
where  p(y;i)  and  x(y}i.)  vary  continuously  with  y.  Since  constant  functions 
are  continuous,  this  result  holds  also  here. 

C.  Exponential  Representations  of  M(x,y) 

The  preceding  definition  of  m(x, y)  is  the  constructive  definition,  the 
one  that  will  at  once  yield  numerical  values  for  the  amplitude  vectors  A(y). 

It  can  be  used  even  when  p  and  x  depend  on  y.  There  is,  however,  a  formal 
definition  of  «(x,y)  that  is  of  great  heuristic  value;  indeed  it  will  lead  us 
to  the  eigenmatrix  representation  of  A(y),  the  central  formula  of  the  present 
study.  This  is  the  definition  that  starts  from  (6.8)  and  visualizes  «(x,y)  as 
given  by  the  same  kind  of  formula  found  in  the  theory  of  the  scalar-valued 
exponential  function.  Thus  let  us  write 

M(x,y)  =  exp[K(y-x)]  (6.13) 

and 

»  ffj(y-x)^ 

exp[K(y-x)]  =  l  - t-j -  (6.14) 

j-0  J’ 

The  operations  in  each  term  of  the  series  (6.14)  are  numerically  possible  and 
convergence  of  the  infinite  series  can  be  established.  Hence  in  principle  the 
exponential  of  the  matrix  K(y-x)  is  computable  arbitrarily  accurately.  It  is 
easy  to  verify  (just  as  in  elementary  calculus)  that 


In  this  case,  then,  exp[K(y-x)]  can  be  evaluated  numerically  quite  readily, 


provided  we  know  the  2qx2q  matrix  E  and  the  2q  numbers  <j,j  =  l,...,2q.  By 
(6.16)  E  and  the  tcj  are  the  eigenstructures  of  K.  That  is,  from  (6.16),  we 
have 


K  E  =  E  <  (6.20) 

which  requires  that  E  =  ( e ^  •••  e^  •••  «2q^  he  thought  of  as  a  matrix 

made  up  of  2qxl  vectors  ej,  j  =  l,...,2q,  each  of  which  satisfies  the 
eigenvector  equation 

K  e.  *  k.  e.  ,  j  =  1 , . . . ,2q  (6.21) 

- J  J  “J 

From  this  we  see  that  is  the  eigenvalue  of  K  associated  with  the 
eigenvector  e^  of  K.  This  fact  about  the  ej  and  tcj  is  central  to  the  present 
study.  We  shall  next  reapproach  (6.20)  from  a  more  physical  direction.  This 
will  allow  us  to  see  the  eigenstructures  and  Kj  as  arising  from  the  local 
reflectance  and  local  transmittance  matrices  comprising  K  in  (6.4). 


‘jrjr 


j  wtj  mvruwv  KUiTU  irvirv  i f\wv *\j wj  irj WV IT JXV NWUilUVl WV ,TTJWW^VriA7V\jWA.V JWWVV 


7.  PHYSICAL  BASIS  OF  THE  EIGENMATRIX  REPRESENTATION  OF  THE  FUNDAMENTAL 
SOLUTION 

We  return  to  the  setting  of  (6.18)  and  (6.20)  to  provide  a  physical  basis 
for  these  formulas.  In  particular  we  may  ask:  what  is  the  physical  basis  for 
the  diagonal  matrix  k  in  (6.19)  and  for  the  vectors  ej  forming  El  Further, 
what  physical  reason  may  be  given  for  the  invertibility  of  E? 


A.  A  Natural  Basis  for  the  Radiance  Amplitude  Vectors 

If  one  plots  the  natural  logarithm  of  N(y;u,v),  for  fixed  u,v,  as  a 
function  of  depth  y,  in  a  deep  homogeneous  natural  water  body,  one  sees  the 
curve  become  essentially  straight  from  some  depth  yQ  downward.  There  is  a 
depth  yQ  for  which  this  is  uniformly  true  for  all  u,v,  u  =  l,...,m; 
v  =  l,...,2n.  Now,  if  the  medium  is  homogeneous  and  infinitely  deep  and  since 
(6.5)  is  a  linear  system,  there  is  the  intuitive  suggestion  that  perhaps  there 

may  be  some  linear  combination  of  the  observable  vectors  A (y)  that  decays  (or 

/ 

grows)  precisely  exponentially  with  depth  y;  and  conversely,  these  purely 
exponentially  decaying  and  exponentially  growing  functions  may  perhaps  be 
linearly  combined  to  yield  A(y).  To  see  where  this  leads,  let  us  postulate 
the  existence  of  2q  distinct  exponential  functions  in  y,  over  the  range 
x  <  y  <  z,  of  the  form 


B  . (y)  =  B  .  (x)  exp[i<  .(y-x)  ] 
J  J  J 

j  =  1 , « • « ,q 


(7.1) 


where  kJ,  j  =  l,...,q  are  2q  distinct  real  numbers.*  They  are  just  as  general 

*  Recall  that  y  in  the  body  of  this  study  is  optical  depth,  so  that  y  =  at , 
where  t  is  the  associated+geometric  depth  and  a  is  the  volume  attenuation 
coefficient.  Hence  the  <"-+are  dimensionless.  They  correspond  to  physical 
attenuation  coefficients  =  cik~. 
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where 


£  *  diag[£+,£  ] 


Furthermore  (7.2)  can  be  written  as 


( 2qx2q  ) 


(7.8) 


jl~(y)  =  B+(y)  <  + 


or  more  compactly  still  as 


(7.9) 


^  B(y)  =  B(y)  £ 


(7.10) 


Let  us  now  return  to  (7.1)  and  observe  that  the  2q  functions  of  y,  B^(y), 
j  *  1 , . • . ,q ,  over  the  depth  range  of  y  in  X[x,z],  are  linearly  independent. 
This  follows  at  once  from  the  fact  that  the  are  pairwise  distinct  and  a 
direct  appeal  to  the  definition  of  linear  independence  of  functions  over  a 
common  domain  (cf.  Courant,  1936,  p.  439).  The  set  of  all  linear  combinations 
of  the  Bj(y)  therefore  forms  a  2q  dimensional  vector  space  of  functions  on 
X[x,z].  We  assume  that  each  radiance  amplitude  A~(y;u),  for  fixed  u,  is  in 
such  a  space.  (This  assumption  will  be  verified  later.)  Then  we  may  write 


X  X 

A~(y;u)  =  £  B+(y)  f+-(u)  +  B.(y)  f.~(u) 

j=l  J  J  j=l  J  J 


(7.11) 


u=l,...,q  ,  x  <  y  <  z , 
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In  this  way  we  have  physically  motivated  the  2q  x  2q  depth  independent  mapping 
F  of  the  exponential  basis  vector  B(y)  into  the  observable  radiance  amplitude 
vector  A(y)  at  each  y,  x  <  y  <  z. 

The  preceding  mapping  can  also  be  postulated  to  go  in  the  reverse 
direction.  Analogously  to  (7.11)  we  may  now  assume  the  existence  of  suitable 
coefficients  ej“(u)  and  ej~(u),  to  be  determined,  so  that 


x  ”1 

B7(y)  =  £  A+(y;u)  et“(u)  +  £  A  (y;u)  e._(u) 

J  u=l  J  u=l  J 


(7.18) 


j  =  1 , .  • . , q  ,  x  <  y  <  z, 


Following  the  pattern  in  (7.12)  and  (7.13)  but  noting  the  difference  in 
running-index  variables  in  each  case,  we  can  write,  with  j  =  l,...,q, 


e+  (u)  =  [et  (u),...,e+“(u)] 

-  1  q 


e  (u)  =  [e  "(u),...,e  “(u)] 

-  i  q 


Then  (7.18)  takes  its  vector  form 

q  q 

B~(y)  =  £  A*(yju)  e+-(u)  +  £  A  (yju)  e“(u) 

u=l  u=l 

In  matrix  form  this  is 


(7.19) 


(7.20) 


i 


Since  A(x)  is  arbitrary,  we  have  the  desired  connection: 

(7.29) 

Thus  we  see  the  linear  combination  coefficients  in  (7.11)  and  (7.18)  are 
connected  directly  with  the  eigenvectors  of  K,  while  the  growth  and  decay 
rates  of  the  B^(y)  are  the  eigenvalues  of  K.  These  are  the  desired 
physical  interpretations  of  E  and  <  in  (7.27). 

We  can  now  work  backwards  to  construct  the  desired  solution  of  (6.5): 

From  knowledge  of  K  we  can  solve  the  algebraic  problem  (7.27)  to  find  the 
matrices  E  and  k.  From  <  we  can  construct  the  B(y)  by  (7.7),  and  from 
F  =  E~l  we  can  construct  the  amplitudes  A(y)  =  B(y)F.  By  (7.28),  (7.29),  and 
(6.18)  we  see  that  indeed  A(y)  is  a  solution  of  (6.5)  with  K  =  E  <  E~ 1 .  Thus 
(7.29)  is  the  fundamental  matrix  solution  yielding  the  desired  radiance 
amplitude  vectors  A(y)  =  A(x)  M(x, y)  at  each  y  in  X[x,z],  and  associated 
initial  amplitude  A(x).  After  some  further  work  on  the  transport  formulation 
of  this  problem,  we  can  translate  the  final  results  of  the  preceding  steps 
into  simple,  elegant  formulas  (cf.  (9.23)  along  with  (16.2)  and  (16.3)). 
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8.  PHYSICAL  FEATURES  OF  THE  EIGENMATRIX  AND  ITS  EIGENVALUES 
The  decomposition  of  the  radiance  field  into  an  upward  (  +  )  and  a  downward 
(-)  set  of  flows  imparts  special  properties  to  the  eigenstructures  E  and  <  of 
the  system  matrix  K  in  (6. A).  By  exploiting  our  physical  image  of  this  two- 
flow  decomposition,  we  can  go  considerably  further  than  standard  differential 
equation  procedures  in  solving  the  system  (6.5). 


A.  Two-Flow  Partition  of  E;  first  form 

Our  first  action  will  be  to  re-partition  the  matricial  set  of  eigenvector 
solutions  e^  in  (6.20): 


E  =  [e  • • •  e  e 

—  —i  — n  — 


-q  — q+1 


associated  with  the  eigenvalues 


(2qx2q) 


(8.1a) 


<lf...,<q,  <q+1,---»<2q 


(8.1b) 


Going  by  the  pattern  (7.21),  let  us  split  each  ej  in  (8.1a)  into  two  qxl 
column  vectors.  Thus  let  us  write 


_+  =  when  j  =  l,...,q 


(8.2) 


=  e^  when  j  =  q+l,...,2q 


(8.3) 


A/Wv  Vw  *v  r.vjvu  »\.  irj  -rfvwvjv^ 


This  is  simply  another  way  of  arranging  the  coefficients  ej_(u)  and 
e^?(u)  in  (7.18).  The  reason  for  this  rearrangement  will  now  become  clear. 


B.  Reversal  Property  of  Eigenvectors  and  Eigenvalues 

Suppose  ej  is  a  2qxl  vector  of  the  form  (8.2)  with  eigenvalue  Kj, 


R 

j  =  l,...,q.  The  reverse  ej  of  ej  is 


defined  as 


U  I 

-q  -q 


i  u 

-q  — q 


(8.4) 


Here  Q  is  the  2qx2q  reversal  matrix.  Next  observe  by  straightforward 
computation  that  Q  has  the  properties 


92  =  I2q 


(8.5) 


q  k  =  -k  q 


(8.6) 


In  view  of  (8.5),  this  last  equation  may  be  written 


9  *  9  = 


(8.7) 


Now  let  e j  be  an  eigenvector  of  K  with  associated  eigenvalue  kj.  Then  by 


def inition 


£  ij  ’  ‘j  -j 


j  =  l » •  •  • » q 


(8.8) 


-  «_  i 

.V. v.  V..-.  .  .  •r.y.  *.  r.  /.V,  -• 


Multiplying  (8.8)  by  Q  on  the  left  and  using  (8.5)  we  find 

2  £  22  £j  =  Kj  2®  j 
which,  by  (8.7)  and  (8.4)  reduces  to 

K  e |  e|  (8.9) 

j  =  l»*»»»q» 

We  conclude  that  if  (ej,  tj)  is  an  eigenpair  of  K,  then  so  also  is 
(e®,  -Kj),  j  ~  l,...,q.  From  this  we  see  that  the  eigenvalues  <j  of  K  come  in 
signed  pairs  kj,  -kj,  j  =  l»...,q.  If  we  re-index  these  eigenvalues  so  that 
in  (8.1b),  Kj+q  =  -lcj  for  j  =  l,...,q,  then  it  follows  that  in  (8.1a), 
ej+(j  =  ej  for  j  =  l,...,q.  Hence  in  (8.2)  and  (8.3)  we  deduce  that,  for 
j  =  1 » • •  •  »q 


++ 

e  . 
“J 


(8.10a) 


-+ 


e . 
-J 


(8.10b) 


In  this  way  we  see  how  the  isotropy  of  the  volume  scattering  function  allows  a 
simplification  of  the  double-direction  superscripts  of  the  eigenvectors. 


ft 


hold  in  such  media.  In  optically  deep  source-free  homogeneous  media  X[x,z] 
with  positive  absorption  coefficient,  a  >  0  (at  some  arbitrary  fixed 
wavelength),  we  expect  that 

(i)  Starting  at  level  x  with  input  A~(x,u)  and  zero  input  A+(z,u),  A  (y»u) 
for  each  fixed  u  eventually  decays  exponentially  with  increasing  depth 
y.  Conversely,  starting  at  level  z,  with  input  A+(zju)  and  zero  input 
A~(x;u) ,  A+(y;u)  for  each  fixed  u  eventually  decays  with  decreasing  depth 

y* 

(ii)  Each  of  the  ♦  and  -  modes  of  decay  in  (i)  has  q  degrees  of  freedom  at  and 

near  the  respective  initial  boundary.  [For  example,  the  sun  can  be 
systematically  raised  above  the  horizon  while  we  are  at  level  x.  For 

each  sun  position,  A “(x)  is  then  seen  to  take  on  a  particular  orientation 

in  its  q  dimensional  vector  space.  Hence  the  full  q-dimensionality  of 
A(y)  for  y  near  x  must  be  employed  to  cover  this  wide  range  of  incident 
radiances. ] 

Now  without  loss  of  generality  we  may  arrange  the  non-negative  members  of 
the  set  of  2q  eigenvalues  ±<j,  j  =  l»«**»q  into  ascending  order: 

<  k2  ^  S  iCq.  From  (i)  we  deduce  that  must  be  positive:  tc x  >  0. 

From  (ii)  we  deduce  that  the  q  eigenvalues  must  be  distinct,  so  that 
0  <  <  k2  <  •••  <  iCq.  Next,  from  linear  algebra  we  know  that  the  pairwise 

distinctness  of  the  2q  eigenvalues  <j  of  K  implies  the  linear  independence  of 
the  set  of  2q  eigenvectors  ej  of  K  (cf.  Franklin,  1968,  p.  73),  so  that  E  has 
an  inverse  E~l  =  F.  Statement  (ii)  actually  implies  a  stronger  property  of  E, 
by  virtue  of  (7.20):  the  set  of  vectors  et,...,e^  and  the  set  ei>...»£q  are 
each  linearly  independent.  Hence  (E+)“l  and  (E-)-1  must  exist  on  physical 
grounds . 


While  the  preceding  properties  (i),  (ii)  and  their  consequences  are  not 
offered  as  mathematical  theorems,  the  reader  may  perhaps  agree  that  for  a  real 


plane-parallel  homogeneous  optical  medium  with  a  >  0  drawn  at  random,  the 
probability  is  zero  that  the  preceding  assertions  about  E  and  K  are  not 
true.  Indeed,  we  can  reverse  the  matter  as  follows.  We  can  say  that  a  choice 
of  a(y,X)  and  o(y;£';5,X)  for  a  homogeneous  medium  is  realistic  provided  p  *  0 
and  provided  the  associated  eigenvalues  -Kj>  j  =  l»...,q  are  such  that 
0  <  Kj  <  ...  <  k  ,  for  all  wavelengths,  all  azimuthal  modes,  and  all  choices 
of  q.  The  resultant  model  is  then  a  realistic  model ,  by  definition.  It  would 
appear  that  a  necessary  and  sufficient  condition  for  a  model  to  be  realistic 
is  that  a(y;X)  and  a(y;5';C,X)  yield  the  inequalities  s(y,X)  >  0  and 
a(y,X)  -  s(y,X)  =  a(y,X)  >  0  for  all  y  and  X. 


E.  Inversion  of  E  by  Partitioning 

It  will  be  useful  to  find  explicit  expressions  for  the  inverse  F  of  E 
defined  in  (7.24).  Our  discussion  in  paragraph  D  above  showed  that  on 
physical  grounds  the  existence  of  E~l  is  practically  certain.  Let  us 
partition  the  2q  x  2q  matrix  F  analogously  to  E  in  (8.11),  so  that 


where  f*  =  [  f  +  (  1 ),...,  f^q )  ] 

— j  j  j 

Then  by  definition,  requiring 


Therefore,  to  find  E-1,  i.e.  to  find  F~ ,  we  may  invert  the  smaller 
matrices  and  their  algebraic  combinations  as  indicated  in  (8.18)  and  (8.19) 
or  in  (8.24)  and  (8.25). 
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*  (y,x)  =  «__(x,y) 


M_+(y,x)  =  M+_(x, y) 

*+_(y,x)  =  *  ,(*.y> 


M__(y,sc)  =  «  (x,y) 


(9.9) 


(9.10) 


(9.11) 


(9.12) 


These  four  statements  constitute  the  interchange  rule  for  M(x,y)  on  X[x,z]. 

It  holds  quite  generally  and  does  not  depend  on  isotropy  or  homogeneity. 
However,  when  homogeneity  is  taken  into  account,  we  can  find  «(z,y)  without 
further  computation  beyond  w(x,y).  The  matrix  Af(z,y)  is  used  for  upward 
evolving  light  fields,  while  #(x,y)  is  used  for  downward  evolving  light 
fields,  as  we  shall  see,  below.  As  for  evaluating  Af(z,y),  suppose  z-y  =  y-x, 
then  Af(z,y)  =  «*J(y,z)  =  if-1(x,y)  =  Af(y,x).  The  latter  follows  from  tf(x,y)  by 
the  interchange  rule. 


C.  The  Downward  and  Upward  Evolving  Radiance  Amplitudes 

From  (6.10)  and  (9.1)  we  can  write  the  mapping  rule  as 


A  (y)  =  A  (x)  M++(x,y)  ♦  A  (x)  M_+(x,y) 
A  (y)  =  A+(x)  h+_(x,y)  ♦  A  (x)  (x,y) 


(9.13) 


(9.14) 


This  is  for  the  downward  evolving  light  field  starting  from  level  x  in  X(x,z], 
x  <  z.  This  may  be  reduced  to  the  numerical  level  by  explicitly  writing  out 
the  components  of  M  (x,y),...,Af _ (x,y).  Thus  from  (9.3)— (9.6)  and  (8.13)  we 
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10.  GLOBAL  INTERACTION  PRINCIPLES 

The  global  interaction  principles  allow  us  to  correctly  include  boundary 
conditions  into  the  fundamental  solution  procedure. 

By  following  the  developments  in  §6c  of  Mobley  and  Preisendorf er  (1988), 
or  in  §7.4  of  H.O.,  Vol.  IV,  we  may  derive  from  the  fundamental  solution  two 
sets  of  global  interaction  principles  for  the  radiance  amplitudes  in  X[x,z]. 
First  the  downward  evolution  set:  for  a  subslab  X[x,y]  of  X[x,z],  x  <  y  <  z, 
we  have 

A+(x)  -  A+(y )  T(y,x)  +  A  (x)  R(x,y)  (10.1) 

A  (y)  =  A+(y)  R(y,x)  +  A  (x)  T(x,y)  (10.2) 

These  statements,  while  written  explicitly  for  inside  the  water  body  part 
X[x,z]  of  X[a,bJ  =  X[a,x]  U  X[x,z]  U  X{z,b],  actually  can  be  phrased  for  any 
subslab  of  X[a,b]  by  replacing  x  and  y  by  other  depth  variables  in  the  range 
[a,b].  For  subslabs  X[x,y]  the  water  body  itself,  we  can  evaluate  the  R  and  T 
matrices  as  follows,  using  the  fundamental  matrix: 

R(y,x)  =  M"j(x,y)  M+_(x,y)  (10.3) 

T(x,y)  =  H _ (x,y)  -  M+_(x,y)  M~|(x,y)  M+_(x,y) 


(10.4) 
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T(y,x)  =  M-'(x,y) 


R(x,y)  =  -«_+(x,y) 


(10.5) 


(10.6) 


Then  the  upward  evolution  set  for  a  subslab  X[y,z]  of  X[x,z],  x  <  y  <  z,  is 


A+(y)  =  A+(  z  )  T  (  z ,  y  )  +  A  (y)  R(y,z) 


A  (z)  =  A  (z)  R(z,y)  +  A  (y)  T(y,z) 


(10.7) 


(10.8) 


Once  again,  these  statements  can  be  extended  to  any  subslab  of  X[a,b]  by 
suitable  replacement  of  y,z  in  (10.7)  and  (10.8)  by  other  depth  variables.  In 
the  case  of  the  water  body  X[x,z]  itself,  the  R  and  T  matrices  are  given  by 


R  (  y ,  z )  =  «-‘(z,y)  M  +(z,y) 


T(z,y)  =  #++(z,y)  -w+_(z,y)  rHz,y)  w_+(z,y) 


(10.9) 


(10.10) 


T(y,z)  =  «-l(z,y) 


R(z,y)  =  -M+_( z ,y)  M"Kz,y) 


(10.11) 


(10.12) 


As  a  consequence  of  the  homogeneity  and  isotropy  of  the  medium  X(x,z],  we 
find  from  the  interchange  rule  and  the  above  M-representat ions  of  R  and  T  that 
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11.  THE  R-INFINITY  FORMULA  AND  SOME  APPLICATIONS 
The  reflectance  matrix  for  an  infinitely  deep  homogeneous  layer  X[x,=] 
plays  a  central  role  in  the  present  study  of  the  light  field.  We  shall  now 
derive  an  expression  for  this  reflectance  and  draw  some  conclusions  about  the 
light  field. 

Now,  from  (10.6)  we  have  an  expression  for  R(x,y),  the  reflectance  of  the 
finitely  deep  layer  X[x,y]  in  X[x,®].  On  letting  y-**=,  using  (9.3),  (9.4)  and 
recalling  that  <!>0,  we  have 


R^  =  lim  R(x,y)  =  lim  -M_+(x,y)  M~|(x,y)  =  ~E  (£  )“ 1 

y-*a>  y-¥oo 


and  in  like  manner,  from  (10.3),  (9.3),  and  (9.5),  we  find 

lim  R(y,x)  =  lim  M~j(x,y)  M+_(x,y)  =  (F*)~lF 
y-*”  y-*<*>  _  R 

— aj 

where  the  last  equality  comes  from  (8.17).  Hence 


(11.1) 


(11.2) 


R  =  ~E  (E+)- 1  =  (F+)“l  F~ 


(qxq  ) 


(11.3) 


The  physical  significance  of  R^,  E±  and  comes  out  on  rearranging 
(11.3)  into  the  forms 


=  -R  E 


(11.4) 


(11.5) 


In  vector  form  these  read 


•  /•-* *  ■  v  v  v  */  v  v  v  v  vV v  v  v  vN,-  v  vV  vV  v  v  v^.'V v 

J’/V/  v v.\\v; v,v.V y >/^v/w>/yltv; y<yywyyyv,v. -v 


4* 


e  .  =  -R  e  . 
“J  -00  ~J 


f.  =  f.  R 
-J  “J  — 


(11.4a) 


(11.5a) 


for  j  *  l,...,q.  Thus  we  see  that  R^  maps  the  "disembodied"  flow  F*  into  F~ , 
and  also  E *  into  -E~ .  Of  course  the  main  interpretation  of  R,,,  is  obtained 
from  (10.7)  on  letting  z-*®.  It  is  clear  by  homogeneity  and  from  (10.5)  that 
as  z->®,  T(z,y)  ■*  0  ,  and  so  in  X[x,®] ,  we  have 


A  (y)  =  A  (y) 
x  <  y  <  ® 


If  we  explicitly  identify  the  entries  of  R^  via 


R  (1,1)  •••  R  (l,q) 


R  (q , 1 )  • * •  R  (q,q) 


then  (11.4a)  and  (11.5a)  state 

q 

e. (r)  =  -  l  Ra)(r,u)  e.(u)  ,  r  =  l,...,q 

J  u=l  J 

q 

f. (u)  =  £  ft(r)  Raj(r,u)  ,  u  =  l,...,q 

J  r=l  J 

where  we  have  used  the  component  relation  defined  in  (8.13)  and  (8.14). 
Another  useful  set  of  relations  comes  from  opening  up  (7.27): 


(11.6) 


(11.7; 


(11.8: 


(11. 9: 
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12.  SURFACE  BOUNDARY  CONDITION:  AUGMENTED  MATRIX  FORM 


The  preceding  discussion  has  shown  that  the  solution  (9.23)  of  (6.5) 
literally  cannot  leave  the  air-water  surface  at  level  x  unless  we  know  the 
initial  amplitude  vector  A(x)  there.  Our  next  main  goal  is  to  compute  A(x). 
This  requires  attention  to  be  redirected  toward  the  air-water  surface  boundary 
conditions  (5.18)  and  (5.19).  These  conditions  show  that  the  mode 
amplitudes  Ap(a;l)  and  Ap(x;i)  just  above  and  below  the  surface  are  coupled  to 
all  other  modes  by  virtue  of  their  interaction  with  the  directionally 
anisotropic  surface.  Hence,  to  proceed,  we  must  now  reinstate  the  presence  of 
the  p  and  l  indexes  in  the  notation.  Since  we  must  consider  all  the  amplitude 
nodes  simultaneously,  we  shall  form  a  vector  from  them,  as  shown  below. 

Now  recall  that  the  number  q  of  components  of  the  lxq  vector  Ap(y;£) 
depends  on  p  and  i,  (cf.  (5 . 10 )-(5 . 15 ) ) :  q  is  either  m-1  or  m,  as  the  case  may 
be.  For  the  present  boundary  condition  calculations  we  can  treat  all  these 
special  cases  in  a  uniform  manner  by  defining  an  augmented  amplitude  vector 
Ap(y;0  of  m  components  for  each  8.  =  0,...,n,  regardless  of  whether  p  is  1  or 
2: 

Ap(y;a)  =  [Ap(y;l;0, . . .  ,A~(y;m;a)  ]  (12.1) 

p  =  1,2;  l  =  0,...,n 

The  mc^  components  of  these  augmented  amplitude  vectors  are  either  zero  or 
N~(y;m,-),  in  accordance  with  (5.6)  and  (5.7). 

We  next  collect  these  n+1  augmented  m-component  vectors  into  one  grand 
m(n+l)  component  vector  for  each  flow: 
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>.m  *  4  .ft*. 


A~(y)  =  (A±(y;0),...,A*(y;n)] 
— p  — p  “P 


P  =  1,2. 


(12.2) 


The  mxm  matrices  tp(a,x;k| l) ,  are  now  gathered  up  into  one  grand 
m(n+l)  x  m(n+l)  matrix  of  the  form  tp(a,x),  where  we  write 


t  (a,x)  = 
“P 


t^(a,x;0|0) 


tp(a,x;0 | n) 


t  (a,x;n|0)  •••  t  (a,x;n|n) 

— P  — p 


(12.3) 


where  p  =  1  or  2.  The  remaining  three  m(n+l)  x  m(n+l)  surface  transfer 
matrices  are  constructed  similarly.  With  these  constructions  (5.18)  and 
(5.19)  become 


A+(a)  =  A+(x)  t  (x,a)  +  A  (a)  r  (a,x) 

— p  — p  — p  — p  — p 

A  (x)  =  A+(x)  r  (x,a)  ♦  A  (a)  t  (a,>-) 

-p  -p  -p  -p  -p 


(12.4) 


(12.5) 
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13.  BOTTOM  BOUNDARY  CONDITION;  AUGMENTED  MATRIX  FORM 


Using  the  augmented  m-component  amplitude  vectors  A“(y;i)  defined  in 


(12.1),  we  can  reformulate  bottom  boundary  condition  (5.22)  in  the  form 


A  (z)  =  A  (z)  r  (z,b) 
“P  “P  “P 


(13.1) 


P  *  1,2 


where  A~(z)  is  now  1  x  m(n+l)  and  fp(z,b)  is  a  m(n+l)  x  m(n+l)  block  diagonal 


matrix  of  the  form 


r  (z,b)  = 
~P 


f  (z,b;0) 
~P 


0  r  (z,b;n) 

— m  — p 


(13.2) 


for  p  -  1,2. 
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14.  IMBED  RULES 


The  imbed  rules  give  the  operators  that  yield  the  amplitudes  at  some 
internal  level  y  of  a  layer  X[x,z]  knowing  the  R  and  T  matrices  for  the  two 
sublayers  X[x,y]  and  X[y,z]  above  and  below  the  level  y,  x  <  y  <  z  (cf.  H.O., 
Vol.  II,  p.  297).  In  the  present  application  of  the  imbed  rule  we  are 
interested  in  finding  the  amplitudes  A~(x)  just  below  the  air-water  surface  in 
X(a,b],  a<x<y<z<b,  where  X[a,x]  is  the  upper  surface  boundary,  X[x,y] 
is  the  water  body,  and  X[z,b]  is  the  lower  boundary  (a  matte  surface  or  a  half 
space).  Given  the  incident  amplitude  A~(a)  =  [A~(a;0), . . . ,A~(a;n)] ,  the 
required  amplitudes  A~(x)  =  [A-(x}0), . .. ,A“(x;n)]  are  given  by 


A  (x)  =  A  (a)  T  (a,x,b) 

~P  ~P  -P 

1  x  m(n+l) 

(14.1) 

A+(x)  =  A  (x)  R  (x,b) 

-p  -p  -p 

1  x  m(n+l) 

(14.2) 

=  A  (a)  R  (a,x,b) 

-p  p 

where  the  complete  transmittance  and  complete  reflectance  operators  are  given 
by 


T 

~P 


R 

-P 


(a,x,b) 

(a,x,b) 


=  T  (a,x)  [I  -  R  (x,b)  R  (x,a)]-1 
— p  —  — P  —  P 

=  T  (a,x,b)  R  (x,b) 

-p  -p 


(14.3) 

(14.4) 


These  follow  from  the  boundary  condition  (12.5)  and  the  global  interaction 
principles  of  §10  written  for  X[a,b]  =  X[a,x]  U  X[x,b].  In  particular  Tp(a,x) 
takes  the  form  tp(a,x),  and  Rp(x,a)  takes  the  form  fp(x,a),  both  of  which  are 
the  m(n+l)  x  m(n+l)  matrices  defined  in  §13,  above.  The  matrix  Rp(x,b)  is 
discussed  in  §15,  below.  Observe  that  by  the  isotropy  of  X[x,z]  and  X[z,b], 
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15.  UNION  RULES 

The  union  rules  give  the  R  and  T  matrices  of  the  union  X(x,b)  of  two 
layers  X(x,z)  and  X(z,b)  (cf.  H.O.,  Vol.  IV,  p.  30),  knowing  the  R  and  T 
matrices  of  the  two  layers.  In  the  present  case  X(x,z)  is  the  water  body  and 
X(z,b)  is  its  lower  boundary.  The  required  rules  follow  at  once  from  bottom 
boundary  condition  (13.1)  and  the  global  interaction  principles  of  §10  written 
for  X(x,b]  =  X[x,z]  U  X[z,b].  They  are 


R  (x,b)  =  R  (x,z)  +  R  (x,z,b)  T  (z,x) 

-p  -p  -p  -p 


R  (x,z,b)  =  T  (x,z)  [I  -  R  (z,b)  R  (z,x)]-1  R  (z,b) 
-p  -p  -p  -p  -p 


(15.1) 


(15.2) 


In  particular  Rp(z,b)  takes  the  form  of  the  m(n+l)  x  m(n+l)  augmented  block 
diagonal  matrix  rp(z,b)  defined  in  (13.2),  while  the  four  matrices  Rp(x,z), 
T  (z,x)  and  R  (z,x),  T  (x,z)  are  m(n+l)  x  m(n+l)  block  diagonal  matrices 

*  r  r 

(augmented  from  (m-1)  x  (m-1)  form  by  adding  zeros  in  the  m*"^1  row  and  m11*1 
columns,  if  necessary)  made  up  of  the  l-mode  matrices  associated  with  the 
water  body.  For  example  we  have 


R  (x,z)  = 

-p 


R  (x,z ;0) 

-p 


R  (x,z;n) 
“P 


(15.3) 


The  mxm  matrix  Rp(x,z;£),  l  =  0,...,n,  is  given  by  (10.9),  and,  as  noted,  is 
augmented  to  mxm,  if  necessary,  for  use  in  (15.3)  with  the  augmented  lxm 
amplitude  vectors  in  (12.1).  The  remaining  three  augmented  matrices  of  the 
water  body  are  assembled  into  block  diagonal  form  similarly. 
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Another  application  of  the  union  rule,  this  time  to  the  union  of  X[a,x] 
(the  upper  surface)  and  X[x,b]  (the  water  body  plus  the  lower  boundary)  yields 


the  matrix  needed  to  find  the  upward  radiance  amplitudes  A*(a)  emerging  from 
the  air-water  surface  of  the  hydrosol.  The  required  rules  follow  from 
boundary  conditions  (12.4)  and  (12.5)  and  the  global  interaction  principles  of 
§10  written  for  X[a,b]  =  X[a,x]  U  X[x,b].  The  resultant  union  rules  are 


R  (a,b)  =  R  (a,x)  + 
“P  ~P 

R  (a,x,b)  =  T  (a,x) 
-P  -p 


R  (a,x,b)  T  (x,a) 

-p  -p 


[ I  -  R  (x,b)  R  (x,a) ]_1 

-p  -p 


R  (x,b) 
“P 


(15.4) 

(15.5) 


Here  R  (a,x),  T  (x,a)  and  R  (x,a),  T  (a.x)  are  the  four  m(n+l)  x  m(n+l) 

—  O  — p  7  — P  '  — p  7 

transfer  matrices  of  the  upper  surface  occurring  in  (12.4)  and  (12.5),  while 
Rp(x,b)  is  the  matrix  found  in  (15.1). 

The  required  upward  emergent  radiance  amplitude  A*(a)  leaving  X[a,b]  are 
given  by 


A  (a)  =  A  (a)  R  (a,b)  1  x  m(n+l)  (15.6) 

-p  -p  -p 
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16.  SOLUTION  FOR  A  FINITELY  DEEP  MEDIUM 


Having  determined  A"(x)  via  (14.1)  and  (14.2),  we  may  now  return  to 
(9.23)  and  find  numerical  values  Of  A~(y;u)  (with  p  and  1  understood)  at  all 
depths  y  in  the  homogeneous  water  body  part  X[x,z]  of  the  complete  medium 
X[a,b]  =  X[a,x]  U  X[x,z]  U  X(z,b].  Indeed,  we  may  now  explicitly  evaluate  the 
initial  amplitudes  a^(x)  in  terms  of  the  two  basic  qxq  reflectance  matrices 
R(x,b)  and  R^  of  the  medium  X[a,b].  First  observe  that  from  (11.4)  we  have 


e.  =  -R  e.  ,  j  =  l,...,q 
J  "J 


(16.1) 


This,  with  the  mode-uncoupled  form  of  (14.2),  allows  us  to  rewrite  (9.19)  as 


a+(x) 

J 


A+(x)  et  +  A  (x)  e . 
“J  “J 


=  A  (x)  R(x,b)  et  -  A  (x)  R  et 
-J  —  -J 


=  A  (x)  [ R( x , b )  - 


(16.2) 


Moreover,  (9.20)  becomes 


a.(x)  =  -A  (x)  R(x,b)  R  e+  +  A  (x)  e. 

J  -  —  — m  — J  —  J 


=  A  (x)  [I  -  R(x,b)  R  ]e . 
-  -  -  -J 


(16.3) 


where  r  and  i  are  understood  in  (16.2)  and  (16.3).  Hence  if  the  optical 
properties  of  the  medium  are  known,  A~(x)  (=  A^(x)  obtained  via  (14.1))  will 
be  the  only  additional  piece  of  information  needed  for  a  full  solution. 


18.  THE  ASYMPTOTIC  RADIANCE  DISTRIBUTION 


When  y  in  (17.3)  becomes  large,  Che  associated  directional  distribution 
of  the  zero  mode  cosine  radiance  amplitudes  takes  a  well-defined  form,  that  of 
the  so-called  asymptotic  radiance  distribution.  The  basis  for  this  is  the 
following.  Recall  that  we  have  arranged  the  distinct,  positive  eigenvalues  of 
each  mode  l  in  ascending  order:  0  <  k^A)  <  <  •••  <  Xq(A). 

Numerical  experiments  with  realistic  models  (cf.  §8D)  invariably  yield 
the  inequalities: 


tc^O)  <  KjCt)  for  all  A=l,...,n. 


(18.1) 


Physically,  this  is  interpreted  as  showing  that  locally  high  curvature  and 
asymmetry  of  the  radiance  distribution  (thought  of  as  a  surface  in  three- 
dimensional  space)  tends  to  be  smoothed  and  reduced  as  depth  y  increases.  An 
intuitive  theoretical  proof  of  (18.1)  can  be  given  along  the  lines  developed 
in  Preisendorf er  (1959).  Now  let  us  multiply  each  side  of  (17.3)  for  the  case 
p  =  1,  A  =  0,  by  eICl^^y_x^  and  then  take  the  limit: 


A  (°°;u)  =  lim  A~(y;u;0)  y  =  a.(x;0)  f*(u;0) 


We  call  the  set  of  numbers  {A~(°°;  1 ) , . . .  ,A“(=>;m) }  or  any  scalar  multiple 
of  this  set  the  asymptotic  radiance  distribution  (of  order  m).  Thus  A~(«>;u) 
is  defined  via  the  zero  mode  (A=0)  cosine  radiance  amplitude  (p  =  1)  of  the 
radiance  field.  All  other  modes  of  the  physical  radiance  distribution  decay 
at  a  greater  depth  rate  than  k^O),  by  (18.1),  and  are  lost  on  the  way  to 
y  =  <=.  Hence  the  asymptotic  radiance  distribution  is  symmetrical  about  the 
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vertical  direction 


.  The  u-dependence  of  f“(u;Q)  defines  the  zenith  to  nadir 
shape  of  the  asymptotic  radiance  distribution.  As  we  have  seen,  fjCujO)  is 
determined  solely  by  the  mxm  system  matrix  K(0)  (for  the  zeroc^  azimuthal 
mode),  which  in  turn  is  defined  by  the  shape  of  the  volume  scattering 
function.  Observe  that  all  directional  information  about  the  initial  radiance 
distribution  via  Aj(x;0)  has  been  lost  in  the  formation  of  a7(x>0)*  Hence  the 
asymptotic  radiance  distribution  is  an  inherent  optical  property  of  an 
infinitely  deep  homogeneous  medium. 
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From  these  one  finds  the  volume  absorption  coefficient  (cf.  H.O.,  Vol.  V, 
p.  247): 


k  (1-R  ) 

CD  03 

D  +  R  D 

—  CD  + 


(19.5) 


which  can  serve  as  a  check  on  the  computations.  Recall  that  the  whole  problem 
began  with  a  and  o  given,  so  that  a  *  a-s,  derived  from  the  initially  given 
data,  can  be  checked  against  a  in  (19.5),  the  end  of  a  long  chain  of 
arithmetic  operations. 

Interestingly ,  from  (19.1)-(19.4)  we  see  that  K(y,±),  R(y,±),  and  D(y,±), 
which  are  apparent  optical  properties  for  small  depths  y,  attain  the  status  of 
inherent  optical  properties  in  the  limit  of  arbitrarily  great  depth. 
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whence 


w(Ay)  =  A~lA2 


(20.7) 


We  are  next  Led  to  find  the  eigenvector  matrix  E  and  eigenvalue  matrix 
k  =  diagt*! , . . . »Kq>  i » • • • ] »  as  follows. 

Evidently  by  inspection  of  (20.2),  the  eigenvectors  of  M(y-x)  are  already 
E±  =  [e^,...,e~],  the  eigenvectors  of  K.  Next,  suppose  we  order  and  then 
label  the  2q  eigenvalues  of  tf(Ay)  as  follows: 


Y„  <  Yn-1  < 

q  q-1 


<  y2  <  Yl  <  1  <  Yl  <  Y2  < 


<  Y+  <  Y+  (20.8) 
q~l  q 


By  (20.2)  we  expect  that 


Yj  *  exp[±KjAy] 


Therefore,  we  write 


ic  .  =  ±( Ay)- 1  in  y~ 
J  J 


,  j  -  1 , . . . ,  q 


(20.9) 


j  -  i»*-»»q 


(20.10) 


Of  course  when  working  with  actual  data  we  will  find  that  we  don't  exactly 
have  YjYj  =  1>  with  the  result  that  we  will  not  exactly  have  <j  =  -Kj.  We 
will  then  adjust  the  found  above  so  that  we  do  obtain  with  this  symmetry 
property.  We  accordingly  set 
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K  .  =  4(tC  .  ~  1C  .  ) 

J  J  J 


j  =  !■>•••  »<i 


(20.11) 


The  E  are  re-computed  for  these  symmetrized  Kj.  Once  E  and  the  < j , 
j  -  l,...,q  have  been  so  found,  we  can  then  go  on  to  compute 


K  =  E  k  E~l 


(2q  x  2q ) 


(20.12) 


which  yields  the  two  estimated  qxq  matrices  p,  t,  discussed  above,  for  each 


mode  l  ■  0,...,n;  whence  a  and  o. 


As  a  check  on  the  preceding  procedure  we  can  construct  the  amplitude 


vectors  A(yj),  j  =  l,...,2q  from  the  estimated  K  in  (20.12)  and  compare  wit 


the  observed  values.  Note  that  the  eigenstructure  < j  in  (20.10)  and  the  E~  of 


Af(Ay)  in  (20.7)  (before  symmetrizing  in  (20.11))  will  reproduce  the  observed 


radiance  amplitudes  A(y^)  in  (20.5)  exactly.  Once  the  symmetrized  and 


their  associated  matrices  have  been  found,  we  will  have  in  effect  estimated 


the  system  matrix  K  of  the  homogeneous  medium,  i.e.,  we  will  know  the  inherent 


optical  properties  a  and  o  of  the  medium.  Let  K  be  the  estimate.  When  new 


incident  light  fields  A'(x)  come  along,  the  associated  new  A(y),  say  A'(y),  at 


any  depth  y  (and  hence  N(y;u,v)  at  the  depth  y),  can  be  obtained  from  the 


mapping  property  A'(y)  =  A'(x)  exp[£(y-x)] 
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APPENDIX  A 


EIGENMATRIX  THEORY  OF  THE  TWO-FLOW  IRRADIANCE  MODEL 

1 .  Introduction 

The  eigenmatrix  theory  of  §7  has  a  simpler  counterpart  in  the  form  of  the 
eigenanalysis  of  the  two-flow  irradiance  model.  We  shall  develop  this  simpler 
version  here,  following  the  main  developments  of  §7  along  a  parallel  track,  as 
far  as  possible.  It  will  be  an  instructive  exercise  on  two  counts!  First, 
the  key  qxq  matrices  and  will  reduce  to  numbers  (because  now  q  =  1)  and 
so  we  will  be  able  to  see  their  physical  constitution  directly;  all  formulas 
can,  if  desired,  be  numerically  evaluated  by  hand,  and  simple  algebraic 
operations  reveal  all  the  inner  workings  of  the  eigenmatrix  theory.  Hence  the 
present  discussion  can  serve  as  an  informative  prerequisite  to  the  main  study 
of  §7.  Second,  the  present  model,  despite  its  simplicity,  is  actually  a  bit 
more  complex  in  the  sense  that  the  local  optical  properties  t+  and  p+,  which 
are  the  present  counterparts  to  the  qxq  matrices  t,  p,  are  in  fact 
anisotropic;  that  is,  unlike  the  radiance  case  of  §7,  we  must  distinguish 
between  absorption  and  backscattering  activities  on  the  upward  and  downward 
flows  on  the  local  level. 

2.  The  Two-Flow  Irradiance  Equations 

The  form  of  the  irradiance  model  we  shall  use  is  that  developed  in 
Preisendorfer  and  Mobley  (1984).  At  any  geometric  depth  y,  x  <  y  <  z, 

*  x-  H(y,±)  =  t  H(y,±)  +  p_H(y,+)  (A2.1) 

ay  -  +• 

where  t+  =  ~(a+  ♦  b+] 
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A2 


and 


a+  =  D+a 


b+  =  D+b 


P+  =  b+ 


We  shall  use  the  p+  and  b+  notation  interchangeably,  according  to  the 
momentary  interpretation  desired:  local  reflectance  or  local  backscatter. 

This  model  has  four  depth-independent  parameters:  the  two  distribution 
factors  D+,  which  describe  the  mean  path  length  of  ascending  (♦)  and 
descending  (-)  photons  through  a  layer  of  unit  thickness;  and  two  inherent 
optical  properties,  the  volume  absorption  coefficient  a  and  the  mean 
backscatter  coefficient  b.  We  shall  work  with  realistic  media  (cf.  §8D), 
i.e..  Chose  for  which  a  >  0,  b  >  0  and  D+  >  D_  >  0. 

The  backscatter  coefficient  b  is  an  inherent  property  in  the  following 
sense.  By  definition,  the  general  depth-dependent  backscatter  function  b(y,±) 
is 


b(y,±)  =  H"i(y,±)  J  d«U)  f  N(y;^')  o(y;c’;0  dfl(s’) 


(A2.2) 


In  the  two-flow  irradiance  model  we  may  adopt  an  average  radiance  over  H+  when 


> 

evaluating  (A2.2).  This  radiance  is  a  form  of  quad-averaged  radiance,  with 

■ft?? 

.  "  w  • 

i 

the  quad  replaced  by 

=+  or  =_. 

Thus  in  (A2.2)  we  may  adopt  radiances  over  5+ 

.hvS> 
,'V' A 

* 

in  the  form 

• 

l 

I 

h(y,>)/2ir 

> 

if  5  is  in  E+ 

i 

N(y,<->  =  < 

h(y,-)/2rt 

(A2.3) 

if  £  is  in  E_ 

:’/v^ 

Under  this  hypothesis 

,  (A2.2)  reduces  to 

• 
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b(y,±)  =  D(y ,±)  b(y) 


(A2.4) 


where  we  have  written 


(A2.5) 


b(y)  s  7-  J  dQ(?)  f  o(y;c';0  dnU') 


or  ^  J  dfl(c)  f  o(y;s';0  daCO 


The  isotropy  of  a  at  depth  y  implies  that  a(y;c'»0  =  o(y for  all  and 
5.  Hence  b(y)  does  not  depend  on  the  direction  of  the  incident  flow  of 
photons  across  the  plane  at  level  y,  and  so  the  alternate  form  of  b(y)  in 
(A2.5)  also  characterizes  b(y).  In  either  of  the  two  formulas,  observe  that 
we  are  finding  a  mean  or  average  magnitude  of  the  backward  scatter  of  photons 
across  the  horizontal  plane  at  level  y,  x  <  y  <  z.  Since  b(y)  in  this  sense 
is  independent  of  the  direction  of  flow  of  the  photons,  it  is  an  inherent 
optical  property  under  the  hypothesis  (A2.3).  Note  that  when  scattering  is 
spherically  symmetric,  i.e.,  when  a(y»C';0  =  s(y)/4iv,  then  b(y)  =  s(y)/2, 
i.e.,  b(y)  is  half  the  volume  total  scattering  function  s(y),  as  expected. 


3.  The  Fundamental  Solution  of  the  Two-Flow  Model 
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We  may  write  (A2.1)  in  matrix  form  as 


^  H(y)  =  H(y)  K 


where  we  have  written 


(A3 . 1 ) 


*■ «.  *  01 

lNi 


A?* 


H(y)  5  [H(y,+) ,  H(y ,-) ] 


(1x2) 


(A3. 2) 


-T  P 

+  + 


“P  T 


(2x2) 


(A3. 3) 


In  this  way  we  see  (A3.1)  as  the  irradiance  counterpart  to  (6.5).  Notice  that 
we  now  have  x  dependent  on  the  upward  (♦)  or  downward  (-)  flow.  Hence  we  do 
not  have  the  local  isotropy  that  is  present  in  (6.4). 

Following  the  development  in  §6,  we  can  write  the  fundamental  solution  of 
(A3.1)  as 


X(x,y)  =  exp[K(y-x)] 


(A3. 4) 


where  K  is  now  the  relatively  simple  2x2  matrix  in  (A3. 3).  It  has  all  the 
properties  of  the  2q  x  2q  matrix  «(x,y)  of  §6.  In  particular,  the  mapping 
property  holds: 


H(y)  =  H(x)  M(x,y) 


(A3. 5) 


We  are  in  effect  working  with  M(x,y)  of  §6  for  the  case  q  =  1.  Keep  in  mind, 
however,  that  the  parallelism  between  the  multimode  theory  and  the  present 
irradiance  theory  is  not  exact,  since  (A2.1)  exhibits  anisotropy  via  the 
property  D+  >  D_  >  0. 
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4.  The  Eigentheory  of  Two-Flow  Irradiance  Fields 

The  purely  exponential  basis  functions  Bj(y)  of  (7.1)  have  their  present 
counterparts  in  two  scalar-valued  functions  B(y,±),  since  we  are  in  effect 
working  in  the  q  =  1  case.  Thus  we  postulate  for  the  irradiance  model  two 
exponentially  varying  functions  B(y,±)  such  that 


B(y,±)  =  B(x,±)  exp(k+(y-x) ] 
x  <  y  <  z 


[  W  •  m~ 2  ] 


(A4.1) 


For  reasons  which  will  become  clear  shortly,  B(y,±)  are  the  eigen-irradiances 
of  the  medium.  Linear  combinations  of  these  are  to  represent  the  observable 
irradiance  fields 


H(y,+)  =  B(y,+)  f++  +  B(y,-)  f_+ 
H(y , - )  =  B( y , + )  f+_  +  B(y ,-)  f__ 


(A4.2) 


where  f++  and  f_+  are  dimensionless  constants.  These  equations  may  be  placed 


in  matrix  form: 


H(y)  =  B(y)  F 


where  we  have  written 


(A4.3) 


B(y)  =  [ B( y , + ) ,  B(y,-) 


f 

-  + 


(1x2) 


(2x2) 


( A4 , 4  ) 


( A4 . 5  ) 
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Conversely,  we  can  represent  the  eigen-irradiances  B(y,±)  as  linear 


combinations  of  the  observable  irradiances  H(y,±): 


B(y,+)  =  H(y,+)  e++  +  H(y,-)  e_+ 
B(y ,-)  =  H(y ,+)  e+_  +  H(y,-)  e 


(A4.6) 


where  e++,  e+_  are  dimensionless  constants.  These  equations  may  be  written 


more  compactly  as 


B(y)  =  H(y )  E 


(A4.7) 


where  we  have  written 


e  e 

++  +- 


e  .  e 


(A4.8) 


Clearly 


F  -  E~ 1 


(A4.9) 


analogous  to  (7.24). 


Corresponding  to  the  law  of  change  of  H(y)  in  (A3.1),  that  for  B(y)  is, 


by  ( A4 . 1 ) , 


B(y)  «  B(y)  k 


(A4.10) 


where 


mm 


VU^Ltfi.k'U'O^. 


k  =  diag(k+,k_] 


(2x2) 


Observe  that  k+  have  units  m-1  because  y  is  now  geometric  depth,  rather  than 
optical  depth  as  used  in  the  body  of  this  report.  Geometric  depth  is  adopted 
here  as  the  more  natural  depth,  since  the  irradiance  model  does  not  have  the 
volume  attenuation  coefficient  a  to  readily  convert  geometric  depth  to  optical 
depth.  (If  it  is  still  desired  to  work  with  y  as  optical  depth  in  (A2.1)  and 
(A4.10),  then  a+,  b+  and  k+  must  be  divided  by  a.) 

Following  the  procedure  leading  to  (7.27)  we  now  find,  via  (A3.1)  and 
(A4.10),  in  the  present  case  that 


K  E  =  E  k  , 


(A4.ll) 


which  is  the  basic  eigenstructure  equation  for  the  irradiance  model. 

Equation  (A4.ll)  contains  two  eigenvector/eigenvalue  statements.  The 
first  may  be  written 


-P  T 


(A4.12) 


In  component  form  this  becomes 


-re  +pe  =  k  e 
+  +•♦  ♦  -+  +  ■*■+ 


e-  +  T-  e-  =  e- 


(A4.13) 


The  two  unknowns  e++  and  e_+  are  determined  by  (A4.13)  up  to  a  common 
factor.  The  first  equation  of  (A4.13)  suggests  that  we  can  set 
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e++  = c 


!<»♦  *  e-+  5  C1(T>  *  k+} 


(A4.14) 


where  Cj  =  1  m.  Thus  e++  has  the  magnitude  of  p  +  and  is  dimensionless.*  The 
second  of  (A4.13),  with  these  values  of  e++  and  e_+  becomes 


-(t+  +  k+)  (x_  -  k+)  +  p+p_  =  0 


(A4.15) 
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We  shall  return  to  this  equation  in  a  moment. 


The  second  eigenvector/eigenvalue  equation  in  (A4.ll)  is 


‘T-  p  + 


-p  x 


(A4.16) 


In  component  form  this  is 


■t  e  +pe  =ke 
+  +-  +  —  -  +- 


-p_  e+_  ♦  t_  e _  =  k_e_ 


(A4.17) 


The  second  of  (A4.17)  suggests  that,  on  using  the  same  C]  as  in  (A4.14), 
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we  set 


ci(t_ 


-  k  ) 


(A4. 18) 
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The  first  of  (A4.17)  yields  (A4.15),  but  now  with  k_  instead  of  k+.  Hence  the 

In  all  subsequent  uses  of  (A4.14)  and  (A4.18),  ct  will  not  explicitly 
appear,  since  its  purpose  is  simply  to  make  the  e++,  e+.  dimensionless.  If 
we  had  adopted  optical  depth  y  in  (A2.1),  the  t+  Ind  p  +  would  be  replaced 
by  the  dimensionless  quantities  x+/ot,  p+/a  and  C[  would  not  be  needed  at 
this  stage. 
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eigenvalues  k+  are  governed  by  the  quadratic  equation 

(k  +  t+)(k  -  x_)  +  p+p_  =  0 


(A4.19) 


or 


k2  +  (t+  -  x_)k  ♦  (p+0_  ~  t+t_)  =  0 


(A4.20) 


From  this  we  see  that  the  roots  k+  of  (A4.19)  satisfy  the  relations 


k+  +  k_  =  t_  -  t+ 


=  [D  -  D_] [a  +  b]  >  0 


(A4.21) 


k+k_  =  p+p_  -  t+t_ 

«  -D  D  a(a  ♦  2b]  <  0 


(A4.22) 


The  roots  themselves  are  given  by 


k+  =  -  T  +  )  ±  [(t+  -  t_)2  -  4(d  +  0_  ~  T  +  T_)]  } 


(A4.23) 


=  -  T+)  ±  [(x+  ♦  T_)2  -  4p+p_]  } 


For  realistic  media  (a  >  0;  see  §8D)  it  follows  from  (A4.22)  that  k+  and 
k_  are  nonzero  and  of  opposite  signs.  From  this  we  conclude  that 


k  <  0  <  k 


(A4.24) 


Under  most  natural  lighting  conditions  (light  from  a  sunny  or  overcast  sky 
entering  a  lake  or  sea),  D+  >  D_,  and  so  from  (A4.21)  we  have  k+  +  k_  >  0;  in 
other  words  |k_|  <  k+.  Therefore,  according  to  this  model,  in  natural  waters 
downwelling  eigen-irradiance  usually  decays  slower  than  upwelling  eigen- 
irradiance. 

Moreover,  if  the  medium  exhibits  back  scattering,  i.e.,  b  >  0,  then  from 
(A4.23)  and  recalling  that  t+  are  negative,  we  observe  that 


6  =  t  +  k  =  x  -  k 
+  +  -  - 


=  4((t+  +  T_)  +  [(t+  +  T_)2  -  4o+0_]^}  <  0 


(A4.25) 


These  inequalities  yield  the  following  physical  interpretations  about  decay 
rates  of  the  photon  streams; 


k  <  a  ♦  b  and  -k  <  a  ♦  b  . 

▼  *r  *r  —  —  *“ 


From  (A4.19)  we  obtain  what  will  turn  out  to  be  key  ratios: 


T+  +  k±  P- 

and,  after  some  rearranging,  we  find  from  this  also  that 


(A4.26) 


(A4.27) 


(p  0  )k  -  62k 

+  —  z  > 

p+p_  -  6  2 


(A4.28) 


The  combination  p  +  p_  -  <5 2  in  (A4.28)  is  never  zero  in  realistic 


backscatter ing  media;  for  we  have  by  (A4.24)  and  (A4.25), 


4  =  p+p_  -  62  =  6(k_-k+)  >  0 


(A4.29) 


Unlike  the  isotropic  setting  of  §7,  the  eigenvalues  k+  of  (A4.ll)  are  not 
of  equal  magnitude;  although,  by  (A4.22),  they  are  of  opposite  sign.  By 
(A4.22)  they  will  be  of  equal  magnitude  and  we  will  have  isotropy  if  t+  =  t_, 
i.e.,  if  in  (A2.1)  we  have  D+  =  D_.  This  condition  unfortunately  is  never 
satisfied  by  irradiance  fields  in  realistic  media;  and  so  we  should  retain 
distinct  values  of  D+  and  D_  in  the  present  irradiance  model.  Occasionally, 
however,  (cf.  Preisendorfer  and  Mobley,  1984,  or  H.O.,  Vol.  V,  p.  64)  it  is  of 
interest  to  consider  the  one-D  case  to  explore  potential  symmetries,  and 
develop  very  simple  light  field  models. 

We  may  summarize  the  preceding  findings  about  the  eigenmatrix  E  in  the 


e++  e- 
_+  e___ 


« 


5  P 


(A4.30) 


The  inverse  F  of  E  may  therefore  be  represented  as 


f+- 


=  A~  1 


p  -6 


-5  p 


(A4.31) 


5.  Eigenmatrix  Representation  of  the  Two-Flow  Model  Fundamental  Solutic 


We  may  now  return  to  (A3. 4)  and  write  it  in  a  form  that  is  parallel  to 


(7.29): 


xf* 


«(x,y)  =  E  exp[k(y-x)]F 


(2x2) 


(A5.1) 


where  E  and  F  are  given  in  (A4.30)  and  (A4.31),  and  k  =  diag[k+,k_],  with  k+ 
as  defined  in  (A4.23).  In  more  detail,  we  have  from  (A5.1)  the  following 
scalar  counterparts  to  (9.3)-(9.6): 


«++U,y) 


M_+(x,y) 


e++  e 


f  +  e  eMy-*>  f 


++  +- 


k. (y-x)  r  k  (y-x) 

e  e  +  J  f  +66  •  y 


\  k^(y-x)  r  k  (y-x)  r 

M+_(x,y)  =  e++  e  +  f+_  +  e+_e  -  7  f 

„  ,  x  k  (y-x)  r  k  (y-x)  £ 

M (x,y)  =  e  e  +  f  +  e e  -  7  f 


(A5.2) 


(A5.3) 


(A5.4) 


(A5.5) 


We  will  occasionally  use  the  e++,...,f _  notation  instead  of  the  P+,6  notation 

since  the  former  notation  acts  as  a  useful  mnemonic.  However,  for  reference, 
we  also  can  write  (A5.2)-(A5.5)  in  the  form 


*++(x>y) 


[p+p_  e 


-  S!  eMy-*>|/a 


y>  -  -  ek-(y-*>]/a 

y>  =  -ptSUfc*<>,-x)  -  J /i 

m _ <xfy)  -  -(«'  ek.(y-x)  -  ot0_  ek-<yx)i/a 


(A5.6) 


(A5.7) 


(A5.8) 


(A5.9) 


It  is  at  once  clear  that 


M++(x,x)  =  M  _(x,x)  =  1 
«+_(x,x)  =  M_+(x,x)  =  0 


(A5.10) 
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Moreover,  M+_(x,y)  and  M_+(x,y)  differ  multiplicatively  only  by  a  constant 
factor  -p+/p_  (=  -b+/b_).  Observe  that  the  interchange  rule  (9.9)-(9.12) 
holds  here  also. 

When  k+(y-x)  is  small,  then  to  first  order  in  k+(y-x)  the  quantities  in 
(A5.6)-(A5.9)  may  be  reduced,  with  the  help  of  (A4.28)  and  (A4.29),  to 


M++(x,y)  =  1  -  t+(y-x)  (A5.ll) 
M_+(x,y)  =  -p_(y-x)  (A5.12) 
M+_(x,y)  =  p+(y-x)  (A5.13) 
H _ (x,y)  =  1  +  t _ (y-x)  (A5.14) 


These  relations  may  be  compared  to  (9.3a)-(9.6a), 


6.  Reflectances  and  Transmittances  for  a  Homogeneous  Lager 

The  formulas  (10.3)-(10.6)  that  convert  the  fundamental  operators 

M++,...,M _  to  reflectances  and  transmittances  R(y,x) , . . . ,T(y,x)  of  a  water 

layer  X[x,y]  can  be  applied  to  the  present  scalar  case  also.  Thus  using 
(A5.6)-(A5.9)  in  those  earlier  formulas,  we  find,  for  x  <  y  <  z, 


at  \  -  £  r  „My-x)  k  C  y— x )  ,  r  k  (y-x)  k  (y-x),_,  .  .. 

R(y,x)  =  -p+6[e  +  J  -  e  -  ][p+o_  e  +  1  -  52  e  -  ]  1  (A6.1) 


■rt  \  -  a  (^J.+k  )(y-x)  ,  k  ( y-x )  .  k  (y-x)., 

Tvx,y)  =  Ae  +  -  1  [p  +  p_  e  +  -52e-/  ]  1 


(A6.2) 


R(x,y)  =  -p_6[eMy_x)  -  ek-(y~x) ] [p+o_  ek-(y’x)  -  ek-(y~x)]-'  (A6.3) 

T( y , x )  =  A ( p+p_  ek>(y_x)  -  62  ek-(y_x)]-!  (A6.4) 
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From  these  relations  we  can  read  off  various  important  physical  properties  of 
the  opticaL  medium.  For  example, 


and 


R(x,y)  *  (p_/p+)  R(y»x)  [=  (D_/D+)  R(y,x)] 

T(x,y)  =  T(y,x)  ■=  e( T-_T. K y-JC ’  T(y,x) 


(A6.5) 


(A6.6) 


The  parenthetical  statement  holds  if  we  adopt  the  assumption  b+  =  D+5,  as  in 
§A2  above.  Thus  we  see  that  upward  and  downward  slab  transfer  properties  are 
equal  if  and  only  if  we  have  lighting  isotropy  (D+  =  D_),  as  observed 
earlier.  When  k+(y-x)  is  small,  (A6.1)-(A6.4)  reduce,  to  first  order  in 
k+(y-x),  to 


R(y,x) 

=  p+(y-*) 

(A6.7) 

T(x,y) 

=  1  +  t_(y-x) 

(A6.8) 

R(x,y) 

=  o_(y~x) 

(A6.9) 

T(y,x) 

=  1  +  t+(y-x) 

(A6.10) 

When  k+(y-x)  is  large,  then 


T(x’y)  =  r4 

T(y , x )  — - 

o+p 


gk_(y-x) 


e-k+(y-x) 


(A6.ll) 


(A6.12) 
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In  the  limit  of  infinitely  deep  media,  we  find,  from  (A6.3),  with  the  help  of 
(A4.27),  (A4.30),  and  (A4.31)  that 


V 

.V 


R_(°°)  =  lim  R(x,y)  =  -6/p  +  =  ~(x++k+)/p+  =  (a++b+  -  k+)/b+ 

y-*as 

=  -p_/ ( x_-k+ )  =  b_/(a_+b_+k+) 


=  -e  e  “l  —  f~ 1  f 
-♦  ++  —  -+ 


The  last  equality,  (A6.15)  (which  stems  from  (A4.7)),  shows  the  formal 
connection  with  the  radiance  case,  namely  (11.3).  The  alternate  formulas 
(involving  a+,  b+,  k+)  are  model  versions  of  exact  relations;  see  H.O., 

Vol.  V,  p.  113.  Furthermore,  from  (A6.1),  again  with  the  help  of  (A4.26), 
(A4.29),  and  (A4.30): 


The  last  equality,  (A6.18),  as  (A6.15),  shows  the  formal  connection  with  the 
isotropic  radiance  case.  Observe  how  (11.12)  and  (11.14)  reduce  to  the 
present  formulas  if  one  momentarily  allows  isotropy  in  the  preceding  equations 
and  relaxes  the  non-commutativity  property  of  matrix  multiplication  in 
section  11  (think  of  the  matrices  as  lxl). 

7.  Solution  for  a  Finitely  Deep  Medium 

We  now  may  assemble  the  pieces  of  the  solution  of  (A2.1)  for  a  light 
field  in  a  finitely  deep  medium  such  as  that  shown  in  Fig.  1.  We  assume  given 
the  set  of  four  irradiance  reflectances  and  transmittances  t(a,x),  t(x,a). 


(A6.13) 

• 

(A6.14) 

(A6.15) 

=  -(x_+k_)/p_  = 

(a_+b_  -  k_)/b_ 

(A6.16) 

m 

*  ~oJ  (  T+-k_)  = 

b+/(a++b++k_) 

(A6.17) 
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r(a,x),  r(x,a)  of  Che  wind-blown  surface  X[a,x],  as  generated  in  Preisendorfer 
and  Mobley  (1985,  1986).  The  optical  properties  of  the  homogeneous  water  body 
X(x,z]  are  specified  by  the  D+,  a,  and  b  parameters  in  (A2.1).  Finally,  the 
reflectance  R(z,b)  of  the  bottom  X(z,b]  is  assumed  specified.  Downward 
irradiance  H(a,~)  is  incident  on  the  upper  surface  Xfa,x]  and  there  are  no 
other  sources  of  flux  on  or  in  X[a,bj  =  X[a,x]  U  X[x,z]  U  X[z,b].  It  is 
required  to  find  H(y,±)  for  all  depths  y,  x  <  y  <  z,  and  also  the  emergent 
irradiance  H(a,+). 

We  start  with  the  mapping  (A3. 5).  Using  (A5.2)-(A5.5)  we  find,  on 
rearranging  the  terms: 


H<y,±>  =  «*<*>  ek.<y'x)  f  +  .  .‘(x)  ek-(y-«>  £_t 


a  (x)  =  H(x,+)  e++  +  H(x,-)  e_+ 
a  (x)  =  H(x,+)  e+_  ♦  H(x,-)  e 


Now  consider  the  composite  medium  X[x,b]  =  X[x,z]  U  X[z,b]  consisting  of  the 
water  body  X[x,z]  and  the  reflecting  lower  boundary  X[z,b].  Let  us  assign  a 
reflectance  R(x,b)  to  X[x,b],  which  we  will  later  show  how  to  evaluate.  For 
the  present,  the  global  interaction  principle  applied  to  X[x,b]  yields  the 
relation 


H(x,+)  =  H(x,~)  R(x,b) 
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The  emergent  flux  H(a,+)  at  LeveL  a  is  given  by  the  global  interaction 
principle  appLied  to  X[a,x]: 


'.-*,sr»v^vT  sr"3^S^'^ 


H(a,+)  =  H(a,~)  r(a,x)  +  H(x,+)  t(x,a) 


(A7.10) 


where  H(x,+)  is  given  by  (A7.3),  and  r(a,x)  and  t(x,a)  are  the  remaining  two 
irradiance  transfer  coefficients  for  the  air-water  surface.  This  completes 
the  solution. 

8.  Solution  for  an  Infinitely  Deep  Medium 

In  (A7.4)  and  (A7.5)  set  R(x,b)  =  R_(<=)  (which  is  the  limit,  as  of 

R(x,b)  in  (A7.6)).  Then  for  x  <  y  <  ®,  (A7.1)  reduces  to 


or  simply 


and 


H(y,±)  =  H(x,-)[1  -  R_(“)R+(“)]  e _ ek-^y  X^f_+ 


H(y,-)  =  H(x,-)Aek-^y 


H(y,  +  )  =  H(y ,-)  R_(o>) 


(A8.1) 


(A8.2) 


(A8.3) 


where  A  is  given  in  (A4.29)  and  R_(”)  in  (A6.13).  H(x, -)  is  given  in  (A7.9) 
with  b  =  °°,  and  H(a,+)  is  given  by  (A7.10).  Thus  we  find 


H(x,-)  =  H(a,-)  t(a,x)  [1  -  R_(<=)  r(x,a)]-1 


H(x,  +  )  =  H(x ,-)  R_(~) 


(A8.4) 


(A8.5) 
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H(a,+)  =  H(a,-)  (r(a,x)  +  t(a,x)  [1  -  R_(»)  r(x,a)]~l  R_(«)  t(x,a)}  (A8.6) 


=  H(a,-)  R(a,=) 


Here  R(a,«)  is  Che  limit,  as  b  ♦  »,  of  R(a,b)  where  by  Che  union  rule  (15.1) 
applied  to  the  union  X(a,b]  of  X[a,x]  and  X[x,b],  we  have 


R(a,b)  =  r(a,x)  ♦  t(a,x)  [1  -  R(x,b)  r(x.a)]-1  R(x,b)  t(x,a)  (A8.7) 
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